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ABSTRACT 

We consider the Helmholtz equation with a discontinuous complex parameter and 
inhomogeneous Dirichlet boundary conditions in a rectangular domain. A variant of the direct 
method of cyclic reduction is employed to facilitate the design of improved multigrid components, 
resulting in the method of CR-MG. We demonstrate the improved convergence properties of this 
method. 


1 INTRODUCTION 


Microwave heating of foods has revolutionised the food processing industry. Effective and 
efficient microwave heating depends very much on a detailed knowledge and understanding of the 
dielectric properties of the food to be processed. This need has given rise to extensive research into 
the dielectric properties of materials; see, for example, Tinga and Nelson [1]. 

Microwave heating can be compared to heating by alternating current. The electric field of 
alternating current changes direction approximately 100 times each second, whereas the microw r ave 
field changes direction approximately 5 billion times each second. The heating effect is 
accomplished by energy transfer to dipoles, most commonly water. Most foods contain between 50 
and 90 percent water. By attempting to follow the very rapidly changing microwave electric field, 
the molecular vibrations of the dipoles are strengthened, thus increasing the temperature of the 
water and hence the food. 

The scalar potential 4> associated with the microwave field satisfies the wave equation 

V V-e/^ = 0, 

which is derived from Maxwell’s equations of electromagnetics. The parameter e describes the 
permittivity of the medium and the parameter the permeability. However, the radiation field in 
a microwave oven varies harmonically in time, and so we look for a solution of equation (1) in the 
form 

<?K X ) i) = e ' ut u ( x ), 

' ^ o 
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Figure 1: A two-dimensional model of a microwave oven. 

where u is a time-independent scalar potential function and u> is the frequency of the microwave 
radiation. By substituting this expression into equation (1), we see that u satisfies the Helmholtz 
equation 

V 2 u + 8 u = 0, 

where 8 := efiu 2 . In general, t and \i are complex numbers, with real parts related to a material’s 
ability to store electrical and magnetic energy respectively, and imaginary parts related to a 
material’s ability to dissipate electrical and magnetic energy respectively. However, the 
permeability of biological materials is close to that of free space, i.e. fi w po = 47T x 10 -7 Hm _1 , 
Hence, since most domestic microwave ovens operate at a frequency of 2450 MHz, we can calculate 
8 for any given permittivity e. 

The oven is represented schematically (in two dimensions) by the rectangular domain depicted 
in Fig. 1. Region 1 corresponds to free space and so e = Co = x 10 -9 Fm' 1 and <5 is a real 
constant in this region. Region 2 corresponds to the heated material and so 8 is a complex constant 
in this region. Energy is fed into the system by a magnetron via the waveguide. Hence, in this 
paper, we consider the solution of the Helmholtz equation with a discontinuous complex parameter 
and inhomogeneous Dirichlet boundary conditions in a rectangular domain. 

We close this section with a plan of the paper. In section 2 we describe the mathematical 
problem and discuss the smoothing abilities of two multigrid smoothers. In section 3 we describe 
the technique of approximate cyclic reduction and show how this can be used to design improved 
multigrid components. Numerical results are presented in section 4 and some concluding remarks 
are made in section 5. 


2 MATHEMATICAL PROBLEM 


Consider the complex two-dimensional Dirichlet boundary value problem 

V 2 u + 8u = 0 in H = Hi U H 2 
s.t. u — g on <9H, 

with data 

8 = 


<5j in subdomain Hi 
82 in subdomain H 2 ’ 


where Hi and H 2 are rectangular subdomains of H (as in Fig. 1). 


(2a) 

(2b) 
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Operator Definitions 


Before attempting to solve this problem by the multigrid method, we need to carefully consider 
the definitions of the discretisation, restriction and prolongation operators. In [2], De Zeeuw 
considers the solution of general linear second order elliptic partial differential equations over 
similar domains. He notes that the rate of convergence of standard multigrid methods often 
deteriorates when the coefficients in the differential equation are discontinuous; he proposes 
matrix-dependent grid transfer operators to overcome these difficulties. However, in our case, the 
discontinuity occurs only in the coefficient of u (viz. 8), and not of Vu. Hence we proceed in the 
following way to define operator V = V(8) in the domain 0, where V can be taken to represent the 
discretisation, restriction or prolongation operator. Firstly, if 8 takes value <5,- in subdomain H, 

(i = 1,2), then we set the value of 8 on the interior boundary between Qi and fl 2 to 
83 := |(<5i + $*)• Secondly, V is defined piecewise by 


V 


V{8\) in 
< 7^(62) in fl 2 . 
V(8z) on df h 


(3) 


In practice, this definition of V, for discontinuous 8 , does not seem to impair the convergence of 
the multigrid algorithm for relevant values of 8. 


Equivalent System of Real Equations 


Consider the discrete analogue of problem (2). Suppose 8 = a + i /3 € <D and g € IR- Using a 
central difference discretisation on a mesh of n x n intervals, the matrix of the discrete system 
Au = f is represented in stencil notation by 



1 


1 

V 

1 


1 


(4) 


where A <E <C (n " 1)2 x C (n-1) \ h := 4 and p := 8h 2 - 4 = (ah 2 - 4) + i/ 3 h 2 . Hence, while most linear 
systems which arise in practice have real coefficient matrices, the discretisation of this problem 
yields a complex linear system. Further applications which give rise to complex linear systems 
include discretisations of the time-dependent Schrodinger equation, inverse scattering problems and 
underwater acoustics. 


A popular approach for solving complex linear systems is to solve the equivalent real linear 
systems for the real and imaginary parts of u. However, the following remarks, due to Freund [3], 
cast doubt on this approach. Suppose that A is a general complex m x m matrix. By taking real 
and imaginary parts, we can rewrite the complex system as the real linear 2m x 2m system 

_ / He A ImA \ ( lieu \ = ( Kef \ 

Bu ~ \XmA —Ke A ) \ -2m u ) [ Xmf J ' 
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It can then be shown that B has eigenspectrum 


<?(B) = { \€<C\\ 2 ea(AA)}, 

which means that a(B ) is symmetric with respect to the real and imaginary axes and hence the 
eigenvalues always embrace the origin. Now if A is complex symmetric (as is the case with (4)), 
then B is a real symmetric matrix with real eigenvalues symmetrically distributed about the origin, 
i.e. B is symmetric indefinite. Therefore the equivalent real system is often harder to solve than the 
original c ompl ex one. 


Smoothing Analysis 


Multigrid smoothing methods are usually basic iterative methods, the properties of which are 
well understood. As the name suggests, the function of a multigrid smoothing method is to reduce 
the rough (high frequency) components of the error as efficiently as possible. This is basically a 
local task and so the smoothing efficiency of a method can be analysed by local Fourier mode 
analysis [4], neglecting interactions with boundaries. The smooth (low frequency) error components 
are reduced on the coarser grids. There is a natural distinction between high and low frequencies 
depending on the type of grid coarsening chosen. Essentially, the low frequencies are those which 
are visible on the coarser grids. In principle, smoothing methods need not be convergent (see [5], 
chpt 7), although in practice most are. 

Consider the discrete analogue of problem (l), An = f, defined on a mesh of n x n intervals. 
Basic iterative methods are based on a matrix splitting A = M — N and are defined by 

Mu (m+,) = Nu^ + f. 


The algebraic error arising from the iterative solution of this system of equations is defined by 
e (m) := u ( m ) _ u and satisfies the equation — Ne^ m \ Denoting the stencils of M and N 

by [M] and [iV] respectively, this equation can be rewritten in stencil form as \M\ e jX +1 ' = [Af] e^. 
Now if we define e^ m+1 ^ := Ae^ m ) and note that the algebraic error can be represented as a 
combination of local Fourier modes 

e(-) = A m e W+*), (0, <f>) € 0 := {(& , 2?) : -s + i < p, q < a}, 

then by substituting this into the stencil representation of the error recurrence we define the error 
amplification factor 

x,n ax - 

— [ M j e ( (,-«+«)• 

The error amplification factor is the factor by which the amplitude of the ( 0 , <f>) Fourier mode is 
multiplied as a result of a single smoothing iteration. Now in the case of standard grid coarsening, 
the sets of smooth and rough frequencies are defined by 

0, :=en(-i,f) 2 , 

0r := 0\0,. 
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Figure 2: Fourier smoothing factors p D for PGS and KACZ. 

Hence the Fourier smoothing factor , which is the worst factor by which all high frequency error 
components are reduced per iteration, is defined by 

p := max |A(0 , <6)1 . 

(«,*)ee r v ' 

Note however that this definition of the smoothing factor is only valid for boundary conditions of 
harmonic type. The influence of Dirichlet boundary conditions can be taken into account 
heuristically (see [6] and [7], for example) in the following way. The error at the boundary is 
always zero and so we define a new set of rough frequencies as 

0^ := 0 r fl {(0 , <f>) € 0 : 0 ^ 0 and/or <j> =/ 0}. 

The corresponding Fourier smoothing factor is defined by 

p D := max |A(0 , <f>) I . 

(0,<J)€0 r D 

This is a mesh-dependent definition. A mesh-independent definition, introduced by Brandt [4], is 
obtained by replacing the discrete set 0 with a continuous analogue, but this is more difficult to 
compute numerically and gives less realistic results in cases where the type of boundary condition 
has much influence. 

There are many possibilities for the choice of smoothing method (see [7], for example), but for 
brevity we consider only two, point Gaufi-Seidel iteration (PGS) and Kaczmarz iteration (KACZ). 
The latter of these two methods, dating back to 1937 [8], is considered here because, when applied 
to the complex linear system Au = f, the method converges for all distributions cr(A) of eigenvalues 
of A. The reason for this is that solving the system Au = f using KACZ is equivalent to solving the 
system AA H v = f with u = A ;/ v (i.e. postconditioning ) using PGS, and the matrix AA H is 
Hermitian positive definite, thus guaranteeing convergence. Applying the smoothing analysis to 
stencil (4), the error amplification factors for PGS and KACZ are 

e '$ 4. e »'^ 

ApCS ~ ~p + e" i0 + e->’ 
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{e ie + e f *)(e" + e*+ + p + p) + 2e , <-®+^ 

Akacz - _ 4 + p p + ( e -^ + e -^)( e -^ + e->> + p + p) + 2e‘( e -' t )’ 

for some p = (ah 2 - 4) + i(3h 2 and (0 , <f>) € 0. Fig. 2 displays contour plots of pg GS and p(? AC z 
plotted as functions of ah 2 and 0h 2 . For fixed values of h and a = Tie 6, as 0 — Xm 8 increases, 

Ppgs increases and Pkacz decreases. Hence we might expect the multigrid convergence rate to 
improve slowly with a KACZ smoother and deteriorate more rapidly with a PGS smoother as /? 
increases . This is borne out in practice. Finally, as a rule of thumb, a good smoothing method has 
a smoothing factor p ® <C In this sense, neither of the two methods considered here is a good 
smoothing method for problem (2). 

3 CYCLIC REDUCTION AND MULTIGRID 

Cyclic reduction (CR) is a direct method of solution for tridiagonal and block-tridiagonal 
systems of linear algebraic equations [9], [10]. For tridiagonal systems which represent 
approximations to 1-D second order ordinary differential equations, CR is as efficient as multigrid 
(MG). For problems in higher dimensions CR becomes too computationally expensive due to fill-in 
within the blocks. However, the design of MG methods in higher dimensions can be facilitated by 
drawing comparisons between MG and CR (see Shaw [11]). 

Approximate Cyclic Reduction 

Consider the system of equations Lu = f. If v is an approximation to the true solution u, then 
we define the error vector as e := u - v and the residual vector as r := f - Lv = Le. Then 
assuming that the error vector e is sufficiently smooth (a condition normally guaranteed by a few 
applications of a smoother in a MG algorithm), the fill-in can be minimised by making accurate 
Taylor expansion approximations of the outlying errors. This method is known as approximate 
cyclic reduction (ACR) [12]. 

Now consider a two-grid method applied to a two-dimensional Toeplitz system. Suppose the two 
grids have mesh sizes h and H — 2h and the fine grid matrix has stencil 

' b 

Lh ~ b a b 

b 

where a and b are scalars. Given an initial approximation v to u, we want to solve the equation 
L h e = r for e to obtain an improved approximation v + e. The method of ACR approaches this 
problem as follows. 

Eliminate the outlying errors in the stencil using neighbouring equations to give 

b 2 

2b 2 0 26 2 

Lh ~ b 2 0 46 2 — a 2 0 b 2 

2 b 2 0 26 2 

b 2 
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This first step of CR has destroyed the band structure of the original five-point operator. Further 
steps of CR would introduce more fill-in, resulting in a relatively inefficient process. Instead, 
assuming the errors are sufficiently smooth, approximate the errors at the NW, NE, SW and SE 
positions (in compass point notation) using accurate Taylor series expansions. This defines the 
ACR-modified coarse grid matrix, which has stencil 


Lh ~ a 


2b 2 

0 

2 b 2 0 Sb 2 — a 2 0 2 b 2 

0 

2b 2 


y 


where a is an arbitrary scaling parameter. From the above information, the definition of restriction 
from the fine grid to the coarse grid can also be gleaned. The ACR-modified restriction operator 
has stencil 

r b i 


Rh 



( 5 ) 


For theoretical considerations it is very convenient to choose restriction and prolongation operators 
which satisfy the relation P 'fa = where R is the adjoint operator of with respect to a 
suitably defined scalar product. However, the adjoint of the five-point restriction operator (5) is not 
a reasonable prolongation (see [13], p. 78). Alternative definitions of the prolongation operator are 
discussed in the following subsection. 


ACR and the Helmholtz Equation 


Consider a two-grid method, with mesh sizes h and H = 2h, applied to the fine grid Helmholtz 
differential operator ChU := V 2 u + 6u. Using a central difference discretisation on a mesh of n x n 
intervals, the fine grid matrix has stencil 


Lh 


_1_ 

h 2 


1 

1 p 1 
1 


where h := £ and p := 6h 2 — 4. Hence a = ^ and b = Now if we choose a — ^ , then the 
ACR-modified coarse grid matrix and restriction operator have stencils 


Lji 


(2hy 


1 0 


l h 


Ry ~ 


1 

0 

4 -^P 2 
0 
1 

1 

1 — p 1 
1 


0 1 
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Therefore the analogous coarse grid Helmholtz differential operator is defined as 
C H u := V 2 u + <5(1 — ^-)u, i.e. ACR suggests solving the Helmholtz equation with a different value 
of 6 on the coarse grid in order to stabilise the MG process. For positive real values of 5 for which 
Lh is indefinite, this corresponds to solving the Helmholtz equation with a smaller value of 6 on the 
coarse grid, thus reducing the indefiniteness of Lh- There are various ways to define the 
prolongation operator. Possibilities include seven-point and nine-point prolongation [14]. However, 
a more effective definition of the prolongation operator for this interface problem is 

4 —4 p 4 

— 4p 3 p 2 — 4p , 

4 — 4p 4 

which is derived from the tensor product of the one-dimensional ACR-modified prolongation 
operator. To extend these ideas to an m-grid process, where h,- is the mesh size of grid H, and 
hi + 1 = 2 h{, we proceed as follows. 



Define := 5 and 6k := <5^_i(l — — ) := 6k - 1 c* (2 < k < m) and pk := 6kh\ — 4. Then the 

differential operator on grid f Ik is defined as 


Cku := V 2 u + 6kU, 


hi 


for 1 < k < m, provided <r = -£■. Therefore, the ACR-modified definitions of the matrix of the 
discrete system on grid f Ik and the restriction and prolongation operators have stencils 


Lk 


R 


fc+i 

k 


1 


i 

i pk i 
i 

i 

i -p k i 

i 


P k i 
P ‘ + ' ~ 


4 —4 pk 4 

-4p fc 3p^ — 4pt 

4 —4 pk 4 


respectively. We call this ACR-modified multigrid process CR-MG. Note that the CR-MG 
restriction operator is similar to the operator naturally suggested by the principle of total reduction 
(see [15] and [16], for example). Further, for Laplace’s equation (i.e. 6 = 0), pt = —4 and the 
CR-MG restriction operator corresponds to half weighting. 


4 NUMERICAL RESULTS 

Consider the complex two-dimensional Dirichlet boundary value problem 

V 2 u -f 6u — 0 in fl = Hi U fU : unit square 
s.t. u = g on dLl, 
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with data 

£ _ f 30 + lOi in Cl 2 : | < x,y < | 

\ 1 in Oi : fl\n 2 ’ 

_ f sin(4y - |)7r on x = 0,§<y<§ 

9 \ 0 elsewhere on dCl 

For convenience, we consider a domain Cl consisting of two concentric squares. The value of 6 in CI 2 
is a typical value calculated from the data in [1], In the following experiment we assess the 
efficiency of the CR-MG algorithm, as described in the preceding section, and compare it with 
standard MG using full weighting restriction and nine-point prolongation. 

The problem is discretised piecewise according to (3) and (4), using central differences on a 
65 x 65 grid. A four-grid method is employed, with standard grid coarsening. This ensures good 
resolution of the inner subdomain CI 2 on the coarsest grid. The multigrid schedule used is the 
V-cycle with two pre-smoothing and two post-smoothing iterations, and LU decomposition with 
partial pivoting is used to solve the defect equation exactly on the coarsest grid. The initial 
estimate is taken to be the zero vector and convergence is measured by log 10 ||r|| 2 , where r is the 
residual vector and ||-|| 2 is the usual Euclidean norm. 

With convergence set to a tolerance of 

l°giolMI 2 < 

the convergence times of MG and CR-MG with PGS and KACZ smoothers were measured and the 
results are displayed in Table 1. All convergence times were measured in seconds on a Sun SPARC 

Table 1: CPU Convergence Times 


time (s) 

PGS 

KACZ 

MG 

22.8 

191.6 | 

CR-MG 

18.5 

155.9 


workstation. We immediately notice that both MG and CR-MG converge much more rapidly with 
a PGS smoother than with a KACZ smoother. This is not unexpected, considering the smoothing 
properties of these two iterative methods. Further, KACZ is a more computationally intensive 
smoother than PGS, having a 13-point stencil as compared to the 5-point stencil of PGS. 

However, most importantly, we find that with both smoothers the rate of convergence of 
CR-MG is significantly faster than that of MG. In fact, with both smoothers CR-MG provides a 
19 percent saving in CPU time over MG. This is a significant saving, especially for larger problems. 
The rates of convergence of MG and CR-MG with a PGS smoother are compared graphically in 
Fig. 3. Both plots are approximately straight lines, a consequence of the grid-independent 
convergence of the multigrid method. 
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MG 

CR-MG 


Figure 3: Convergence of MG and CR-MG with a PGS smoother. 


5 CONCLUDING REMARKS 

In this paper, attention has been focussed on improving the design of the standard multigrid 
method with respect to a particular problem, namely the complex-valued microwave oven problem. 
By drawing a comparison with the direct method of cyclic reduction, improved discretisation, 
restriction and prolongation operators have been designed, resulting in savings of up to 19 percent 
in CPU time used. 

Only two smoothing methods have been considered here, point GauB-Seidel and Kaczmarz. 
However, there are many more robust smoothers, such as alternating damped Jacobi, alternating 
symmetric line GauB-Seidel and incomplete LU decomposition. These methods, and many more, 
have been summarised and analysed in detail in [7], Improvements in the convergence properties of 
the modified multigrid method (CR-MG) will almost certainly be realised by using such smoothers. 

Finally, attention in this paper has been restricted to the microwave oven problem, although the 
ideas presented here can be extended to other problems. For example, in [11], these ideas were 
applied to the convection-difFusion equation and it was shown that approximate cyclic reduction 
can be used to define the ideal quantity of coarse grid artificial viscosity and the direction in which 
it lies. 
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